## 
## Author: Kai Arzheimer
## URL: http://www.kai-arzheimer.com/tag/PoolingPolls
## 
######################################################################
## 
## 
## 
######################################################################
## 
### Change log:
## 
## 
######################################################################
## 
### Code:

library(rjags)
library(coda)
library(foreign)
library(R2jags)

# Read Statafile
series <- read.dta("../data/mergedpollsseries.dta")

# Create vote shares from absolute numbers
befragte <- series$v7
anteil <- series[c("cdu","spd","gruene","fdp","linke")]
anteil["sonstige"] <- 100 -rowSums(anteil)
anteil <- round(anteil / 100 * befragte)
befragte <- rowSums(anteil)

# Extract pollsters
company1 <- series$company1
company2 <- series$company2
company3 <- series$company3
company4 <- series$company4
company5 <- series$company5
company6 <- series$company6

pollster <- as.numeric(series$pollster)

# Extract week
hypweek <- series$hypweek + 1

# Extract number of weeks
J <- length(levels(factor(hypweek)))

# Calculate starting values for random walk

woche0sgfl <- series[c("spd","gruene","fdp","linke")][series$hypweek==0,]
woche0cdu <- series[c("cdu")][series$hypweek==0,]
woche0sgfls <- woche0sgfl
woche0sgfls[,"sonstige"] <- (100  - (rowSums(woche0sgfl) + woche0cdu))

partystartlogits <- matrix(data=rep(NA,12),ncol=2)
for(i in 1:5) 
    {
        # Start in line 2 with min (1st col)
        partystartlogits[i+1,1] <- log(mean(((woche0sgfls -2) / (woche0cdu +2))[,i]))
        partystartlogits[i+1,2] <- log(mean(((woche0sgfls +2) / (woche0cdu -2))[,i]))
    }
    
# Modules
load.module("glm")
load.module("dic")


rpr2 <- jags(model="random-pollsters.bug",data = list("anteil" = anteil, "befragte" = befragte,"hypweek"=hypweek,"J"=J,"pollster"=pollster,"partystartlogits"=partystartlogits),n.chains=3,parameters.to.save=c("T[1:39,1:6]","coalmaj[39]","fdpin[39]","redgreenmaj[39]"))

# Save results

btw <- autojags(rpr2,n.iter=5000)
btwsummary <- btw[[2]][[10]]

save(btw,btwsummary,J,file="../data/btwsummary")

######################################################################

